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Q ! Abstract 

^ \ We have studied precursors of the global failure in some self-organised critical models of 

sand-pile (in BTW and Manna models) and in the random fiber bundle model (RFB). In 

-(-i . both BTW and Manna model, as one adds a small but fixed number of sand grains (heights) 
c^ ■ 

2 I to any central site of the stable pile, the local dynamics starts and continues for an average 

'^ ' relaxation time r and an average number of topplings A spread over a radial distance ^. We 

S . find that these quantities all depend on the average height hav of the pile and they all diverge 

! as hav approaches the critical height he from below: A ~ (he — hav)~^ , r r^ (he — hav)~'^ 

^ ; and i ^ [he- hav)'" ■ Numerically we find 5 ~ 2.0, 7 ~ 1.2 and v ~ 1.0 for both BTW 

P^ ■ and Manna model in two dimensions. In the strained RFB model we find that the break- 

Q ! down susceptibility x (giving the differential increment of the number of broken fibers due 

(^ \ to increase in external load) and the relaxation time r, both diverge as the applied load or 

C^ ■ stress a approaches the network failure threshold ac from below: x ~ (c"c — cr)^^/^ and r ~ 

I . {oc — cr)^^/^. These self-organised dynamical models of failure therefore show some definite 

C , precursors with robust power laws long before the failure point. Such well-characterised 
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^ ' precursors should help predicting the global failure point of the systems in advance. 
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I. Introduction 

In a sand-pile, whenever the local slope at the surface of the pile exceeds the angle of repose, 
avalanches take place and the sand grains move to the neighbouring sites. If the local slope of 
these neighbouring sites increase, in turn, beyond the angle of repose, avalanches continue. 
Otherwise the dynamics stops until another sand grain is added to the pile. The system 
finally attains a self-organised state where extra grains, when added, get out of the system 
through successive avalanches from its boundary. Models of sand-piles have been developed 
to study such self-organisation. Bak, Tang and Wiesenfeld (BTW) ||I], |^ introduced the 
random height sand-pile model where height units are added randomly at any site at a 
constant rate and a site topples when its height equals an integer threshold value hth (= 4 
for square lattice, for example). Whenever any site topples, the local height becomes zero 
there and the height is locally conserved by equal sharing among the nearest neighbours 
(four in number for square lattice) and the neighbours get one unit of height added to theirs'. 
The boundary of the system is completely absorbing. As more and more grains (heights) 
are added slowly to the system, the average height hav of the system gradually increases 
and attains a critical height he (equivalent to the angle of repose of the sand-pile), beyond 
which the growth of average height stops as the further addition of grains at any site causes 
successive avalanches or failures of all sizes. These happen due to the long-range correlations 
developed and the additional grains finally get out of the system through its boundaries. The 
self-organised state here becomes critical as it involves power law behaviour in avalanche 
size distribution and the corresponding lifetime distribution. Extensive numerical checkings 
confirmed this self-organised critical behaviour in both two and three dimensional sand piles 
0, ||. Later, Manna introduced a two state and stochastic version of the BTW model 
where the threshold height has been chosen to be two {hth = 2). The toppling at any site 
reduces the height there to zero and the toppled heights add to the height of any stocastically 
chosen site among the four neighbouring sites of the toppled one. Here also, with constant 
addition of sand grains, the system gradually reaches again a critical state and there the 
avalanche size distribution and the corresponding lifetime distribution again follow similar 
scaling behaviour. However, the exponents for the Manna model seem to be different H, |^ 
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from those of the BTW model. A similar self-organising dynamics is also seen in a strained 
random fiber bundle (RFB) model 0, ||, ||, |10|, [ll|, |l^ where A^ fibers are connected in 



parallel to each other and clamped between their two ends. The strength of the individual 
fibers has a random distribution (white, Gaussian or otherwise). Under a load F, a fraction 
of the fibers fail immediately whose strengths are less than the stress a (= F/N). After this, 
the total load of the bundle redistributes globally as the stress is transferred from broken 
fibers to the remaining unbroken ones. This redistribution causes secondary failures which 
in general causes further failures and so on. After some typical relaxation time r (dependent 
on a), the system ultimately becomes stable if the applied stress a is less than a critical 
value (Tc, beyond which all the fibers break and the network fails completely. Although the 
RFB model is not a self-organised critical one (as the failure state at a > cTc is not critical) , 
it has some self-organising dynamics (stress redistribution for a < o"c) similar to the earlier 
ones and is very simple to tackle analytically. The studies of these self-organising model 
systems and their scaling behaviour have been extremely useful in analysing the statistics of 



fracture and breakdown in real materials, including in earthquakes U, [T^, [TJ . 

An obvious question arises: Are there any precursors or prior indications which can tell 
how far a (slowly) growing sand-pile or a gradually strained fiber bundle is away from its 
global failure point? The study of precursors in self-organised systems was initiated by 
Acharyya and Chakrabarti ||l5l. Here, the global failure is identified as the system spanning 
avalanche occurring at hav = he. They tried to study the response of BTW model to pulsed 
addition of grains (heights) in two and three dimensional sand-piles, where, 'pulse' means a 
fixed number of grains, added at any site to trigger the dynamics locally in time and space. 
Adding a pulse of heights at any site of a stable pile (where toppling had stopped), they 
measured the response of the system in terms of the number of affected or toppled sites (A) 
and the corresponding response or relaxation time (r) at various average heights (hav) of 
the system. They observed that both A and r diverge as hav approaches the critical height 
he- They also estimated the exponents involved in the power laws for these divergences. 
However, these estimates for the exponent values were not quite accurate due to the small 
system sizes considered and strong pulses applied. Similarly, the breakdown susceptibility 



T5| of the RFB model was studied by measuring the increment in the number of broken 



fibers with the increment in the stress a 16 . It was seen that this differential increase in 



the number of broken fibers due to infinitesimal increase in stress a, diverges as the stress a 
approaches the global failure threshold ac. 

In this paper, we have studied several precursors in the models of sand-piles and random 
fiber bundle. We have studied the response of sand-pile models (both BTW and Manna 
model) to pulsed addition of sand grains (heights; for unit time or unit pulse width), where 
the applied pulse strength is negligible, so that the statistical state of the system is not 
perturbed significantly by the applied pulse. We have identified three parameters, namely, 
the total number of topplings (A), the corresponding relaxation time (r) and the correlation 
length (^); all of which diverge as the average height {hav) of the pile approaches the critical 
height {he). The values of the exponents for the variations of these quantities (A, r and 
^) with hav near he have been estimated accurately. In fact, the estimated value of the 
critical height or the location of the catastrophe point h^ extrapolated separately from the 
growing (precursor) values of A, r and ^ (for hav values below he), agree quite well with the 
previous direct numerical estimates [l^ for the same. In the RFB model, we have studied 
the breakdown susceptibility (x) and the response time (r) required for the bundle to become 
stable when an initial load or stress a (< a^) is applied on it. Both x and r diverge as a 
approaches Uc. The growth behaviour of these precursors for a below o"c and the possibility 
of their extrapolations for estimating the failure point ac of the network is discussed. 

II. Precursors in the BTW model 
a) Model 

Let us consider a BTW model on a square lattice of size Lx L. At each lattice site {i,j), 
there is an integer variable hij which represents the height of the sand column at that site. 
A unit of height (one sand grain) is added at a randomly chosen site at each time step and 
the system evolves in discrete time. The dynamics starts as soon as any site {i,j) has got a 
height equal to the threshold value {hth= 4): the site topples, i.e., hij becomes zero there. 



and the heights of the four neighbouring sites increase by one unit 

hij -^ hij — 4, hi±ij -^ hi±ij + 1, hij±i — > hij±i + 1. (1) 

If, due to this toppling at site {i,j), any neighbouring site become unstable (its height reaches 
the threshold value), they in turn follow the same dynamics. The process continues till all 
sites become stable {hij < hth for all («,j)). When toppling occurs at the boundary of the 
lattice (four nearest neighbours are not available) , extra heights get off the lattice and are 
removed from the system. 

With a very slow but steady rate of addition of unit height (sand grain) at random sites 
of the lattice, the avalanches get correlated over longer and longer ranges and the average 
height (hav) of the system grows with time. Gradually the correlation length (^) becomes 
of the order the system size L. Here, on average, the additional height units start leaving 
the system as the system approaches toward a critical average height hc{L) and the average 
height remains stable there (see Fig. 1). Also the system becomes critical here as the 
distributions of the avalanche sizes and the corresponding life times follow robust power 
laws 1^, |]. In fact, a finite size scaling fit hc{L) = /ic(oo) + CL^^^'^ (obtained by setting ^ 
~ I hc{L) — hc{oo) |~^= L), where C is a constant, with z/ ^ 1.0 gives he = hc{oo) ~ 2.124 
(see inset of Fig. 1). Similar finite size scaling fit with u = 1.0 gave hc{oo) ^ 2.124 in earlier 



large scale simulations ||T^ 



b) Simulation studies for pulsed perturbation 

We have taken random height BTW systems on square lattice of different sizes (L = 100, 
200 and 300). At a fixed value of L, for any pile configuration at an average height hav, 
when all sites of the system have become stable (dynamics have stopped), a fixed number of 
height units hp = 4 (pulse of sand grains) is added at any central point of the system. Just 
after this addition, the local dynamics starts and it takes a finite time or iterations to return 
back to the stable state (hij < 4 for all {i,j)) after several toppling events. For each value of 
hav{< he), we take about 10^ initial configurations and this response or relaxation time has 
been noted for each of them. The average relaxation time r is obtained taking averages over 



all configurations and is seen to diverge as hav approaches the critical height he (see Fig. 2 
(a)). Near h^ r follows a power law r ~ {he — hav)~'^ ■, where 7 = 1.2. The plot of t~^/'^ 
with hav is a straight line with negative slope. Extrapolating the straight line and locating 
the vanishing point of r"^/''' one can estimate the critical point he = 2.13 ± .01 (see inset of 
Fig. 2 (a)) which is very close to the previous numerical estimate he = 2.124 fi% . 

Another response parameter, the average size of the damage (A), i.e., the average number 
of topplings (after the addition of pulse) has been measured as follows: the number of 
topplings for each configuration at each value of hav is noted and averaged out over the 
initial configurations (about 10^ in number). Thus the average A for that value of hav is 
estimated and this is also seen to diverges as hav -^ he (see Fig. 2 (b)). Near the critical 
point, we find A ~ (/i^ — hav)^^, where 5 = 2.0. The plot of A^^/'^ versus hav gives a straight 
line with negative slope (see inset of Fig. 2 (b)) which can again be used to estimate he 
(= 2.12 ± .01) after extrapolating the straight line up to the vanishing point of A^^/^. 

We have also measured the correlation length ^ of the system during the same experiment. 
When the pulse is added at any central point {io,Jo) of the system at some hav, toppling starts 
there and gradually it moves toward the boundaries . We have marked the farthest affected 
site {if,if) (where at least one toppling has occurred due to the pulse) with respect to the 
central site (^o? Jo) where the pulse had been added. Clearly, the average (over configurations) 
distance between the central and the farthest affected sites (|(«0) Jo)~(V' J/)l) i^ ^ measure of 
the correlation length of the system at that hav This correlation length ^ is seen to diverge 
as hav -^ he (see Fig. 2 (c)) following a power law C, ~ {he — hav)^" , where z/ = 1.0. The plot 
Qf ^-i/i' versus hav (see inset of Fig. 2 (c)) is a straight line. The vanishing point of ^"^/'^ 
gives an estimate of the critical point he and we find he = 2.13 ± .01. This is also close to 
the previously estimated critical value. 

III. Precursors in the Manna model 
a) Model 

We consider now the Abelian Manna model on a square lattice of size L x L, where the 
sites can be either empty or occupied with unit height i.e., the height variables can have 



binary states hij = 1 or hij = 0. A site is chosen randomly and one height is added at that 
site. If the site is initially empty, it gets occupied: 

hi J -^ hi J + 1, (2) 

If the chosen site is previously occupied then a toppling or 'hard core interaction' rejects 
both the heights from that site: 

and each of these two rejected heights stochastically chooses its host among the 4 neigh- 
bours of the toppled site. The toppling can happen in chains if any chosen neighbour was 
previously occupied and thus cascades are created. After the system attains stable state 
(dynamics stopped), a new site is chosen randomly and unit height is added to it. Thus the 
system evolves in discrete time steps. Here again the boundary is assumed to be completely 
absorbing so that heights can leave the system due to the toppling at the boundary. 

With a slow rate of addition of heights at random sites, initially the average height of 
the system grows with time and soon the system approaches toward a critical average height 
he, where the average height stabilizes and does not change with further addition of heights 
(see Fig. 3). The critical average height he has a finite size dependence and a similar finite 
size scaling fit hc{L) = hc{oo) + CL^^^'^ gives z/ ~ 1.0 and he = hc{oo) — 0.716 (see inset of 



Fig. 3). This is close to an earlier estimate he — 0.71695 [0, made in a somewhat different 
version of the model. The avalanche size distribution has got power laws similar to the BTW 
model, at this self-organised critical state at hav = he. However the exponents seem to be 
different H, |^, compared to those of BTW model, for this stochastic model. 
b) Simulation studies with pulsed perturbation 

We have considered Manna model on square lattice of different sizes {L = 100, 200 and 
300). At a fixed value of L, for any pile configuration at an average height hav, a fixed 
number of heights hp = 2 has been added at any central point of the stable pile (for which 
dynamics had stopped). Just after the addition, the local dynamics starts and it takes 



a finite time (iteration number) to return back to the stable state {hij < 2 for all (i, j)) 
after several toppling events. For each value of hav{< he) this response time for each pile 
configuration has been noted and the average relaxation time r is obtained from the average 
over 10^ different configurations. Near critical point r is seen to diverge (see Fig. 4 (a)) as 
hav approaches the critical height he with a power law r ~ {he — hav)^'' , where 7 = 1.2. 
The plot of T'^^"' with hav is a straight line (see inset of Fig. 4(a)) with negative slope. 
Extrapolating the straight line and locating the vanishing point of r'^^"' we have estimated 
the critical height as he = 0.72 ± .01, which is very close to the previous numerical estimate 
he = 0.716 for this model (see inset of Fig. 3). 

The size of the damage, i.e., the total number of topplings (after the addition of pulse) 
has also been measured for the above cases. The average (over about 10^ configurations) 
number of topplings A also diverges as average height hav approaches the critical height he 
and near critical point A grows as A ~ {he — hav)~^, where 6 = 2.0 (see Fig. 4 (b)). The 
plot of A^^/^ versus hav gives a straight line which can be used to estimate he (= 0.72 ± .01) 
after extrapolation (see inset of Fig. 4 (b)). 

The correlation length (^) of the system has been measured following the same procedure 
as in the BTW model, described in the previous section. The average (over about 10^ 
configurations) correlation length ^ again diverges (see Fig. 4 (c)) as hav -^ he and near 
critical point ^ follows the power law ^ ~ (he — hav)"", where z/ = 1.0 . The plot of ^^^^'^ 
versus hav is a straight line with negative slope and the vanishing value of ^"^/'^ estimates the 
critical density he = 0.72 ± .01 (see inset of Fig. 4 (c)) which is again close to the estimated 
critical density from direct numerical study. 

IV. Precursors in the random fiber bundle model 
a) The model 

We consider a RFB model containing N elastic fibers clamped at two ends, where the 
failure stress of the individual fibers are distributed randomly and uniformly within and 1 
(white distribution). Global load sharing is assumed and the applied load on the bundle is 
democratically shared among the existing intact fibers of the bundle. With the application of 



any small load F (= a N, with a <^ 1) on the bundle, an initial stress a sets in. At the first 
step, aN number of fibers are broken off, leaving Nui{a-) = (1 — o-)N number of unbroken 
fibers. After this, the applied force is redistributed uniformly among remaining intact fibers 
and the stress (per fiber) is then readjusted to a value F/[Nui{a)] = 0-/(1 — a). With 
this new readjusted stress, some extra fibers for which the strengths are below the above 
readjusted stress fail and the total number of broken fibers increases to a value N[(7/{1 — a)], 
leaving Nu2{o-) = [1 — o'/(l — <^)]N unbroken fibers. This in turn readjusts the stress again 
and induces further failure giving rise to a recursive relation: 

Wn(^) = 1 ^, (4) 

for the fraction u of unbroken fibers at the n-th and (n — l)-th iteration for stress a. This 
dynamics of successive failure propagates therefore in (discrete) time until Nun-i{cr) — 
Nun{<y) < 1, or the successive stress readjustments make so little change that even one fiber 
cannot be found in the network having strength between the successive readjusted value. 
For an infinite (A^ -^ oo) fiber bundle, we denote the fraction of unbroken fibers here by the 
fixed point value u*{(t). The critical stress Uc is determined by that a above which there is 
no fixed point and Mn(o') ^ as n — > oo. Because of the above simple recursion relation (4) 
for u, in the uniformly distributed RFB model, we can easily analyse the asymptotic features 

of its dynamics. The differential form of the above recursion relation (4) can be written as 

du {u^ — u + a) 

dn u 

The fixed point value of u is obtained by setting du/dn = 0. This gives 

u^ = I + {a. - cry/\ (6) 

where cTc = 1/4. The other root is neglected here as it is unstable (see equation 7). Expanding 
the equation (5) near the fixed point value (6) of u, we can write u = u* + e, and 

de e{2u* - 1) _ ^^^^_ -^/^ 



dn u* 



K^c-ay^% (7) 



as 0" —i> (Tc, which gives 



where 



Un = u* + const. exp(— n/ro), (8) 



ro = ^(a,-cr)-V2. (9) 

b) Study of the precursors 

We have simulated the RFB model with a very slow but steady increase of initial stress 
cr on a bundle containing N fibers (A^ ~ 10^). Application of some small initial stress a 
(= F/N) triggers the dynamics by breaking off a fraction (1 — u„) of fibers, and global 
readjustment of the stress causes further failures (u„+i < m„). As mentioned before, after a 
few steps or iterations, when N[un--i{cr) — n„(cr)] < 1, the dynamics stops and the bundle 
becomes stable. We note this relaxation time r required for the stabilisation. For each 
(initial) stress a we start afresh with the intact bundle and note the relaxation time for each 
a. The observation continues until we reach the threshold stress dc (= 1/4), above which 
the bundle fails totally (see Fig. 5). The relaxation time r is seen to diverge as a ^ ctc 
following a power law r ~ tq ~ (ac —a)^^^"^ (see inset of Fig. 5) which can be explained 
easily using equation (8). 

Similar studies have been made for the breakdown susceptibility x— dm /da, where m 
= N{l—u*{a)) is the total number of fibers broken finally by stress a (see inset of Fig. 5). 
One finds x ~ (o"c — cr)^^^, in agreement with the previous observations ||ll], |16[. This can 
be easily explained from solution (6). 

V. Summary and concluding remarks 

In all the three dynamical models of failure we have considered here, we find that long be- 
fore the occurrence of global failures, the growing correlations in the dynamics of constituent 
elements manifest themselves as various precursors. The number of topplings A, relaxation 
time r and the correlation length ^, in both BTW and Manna model, grow and diverge 
following power laws as the systems approach their respective critical points he from below: 
A ~ (/ic — hav)~^, T ~ {he — hav)~^ and ^ ^ {he — hav)'" ■ For two dimensional systems, we 



find numerically here 6 ~ 2.0, 7 ~ 1.2 and z/ ~ 1.0 for both BTW and Manna model. We 
could not thus detect any significant diff'erence in the power laws for these precursors. We 
also could not detect any significant finite size effect in these precursors. Though this size 
independence of the quantities we studied look quite unnatural at first sight, there are strong 
reasons. Basically, we study the behaviour for hav < h^ the precursor behaviour, where ^ 
is necessarily finite. As we add here the tiny pulse at some central site of a relatively large 
system, the boundary effect can not be really felt because of the smallness of ^ compared 
to L for most values of hav This explaines the lack of finite size effect in our precursor 
studies (which of course is clearly manifest when we check our model results at hav = hc{L)). 
It may also be noted that since for hav near h^ in our system, ^ becomes of the order of 
L, at hav = hc{L), our results suggest A ~ L^° and r ~ L^'^. This in fact supports the 
earlier analytic result for hav = hc{L) for large but finite systems, as obtained by Dhar fl^ . 



Generally, if we write A ~ ^'^f , we then get dj = 5 /u ~ 2.0 for the fractal dimension of the 
avalanche clusters. 



Apart from the previous attempts ||15|| , an indirect study in the fixed energy sand pile 
(FES) model [|1^] also indicated similar power law behaviour away from the critical point 
(essentially for hav above he). In the BTW-FES model, the observed exponent values for 
r and ^ differ significantly from those of ours'. However for the Manna-FES model, these 
exponent values are close to our estimates. The FES version of the models are somewhat 
different by construction and the discrepancies in case of BTW-FES estimates (compared 
to ours') seem to be physical in their origin. Due to lack of stochasticity, BTW model can 
stabilise in several 'metastable' states (above ho) and non-universality occurs because of 
different initial conditions. This can also be seen from the difference in the estimate of the 
critical point he in the BTW and BTW-FES models; as mentioned before, no such difference 
in the he estimate seems to exist for the Manna and Manna-FES models. This difference 
in the h^. values for the BTW case might explain the difference in the exponent values we 
obtained (for hav < he) and those obtained for the corresponding FES model (for hav > he). 

For the random fiber bundle model, we find that the breakdown susceptibility x (giving 
the increment in the number of broken fibers for an infinitesimal increment of load on the 
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network) and the corresponding relaxation time r (required for the network to stabihse, after 
successive failures of the fibers), both diverge as the external load or stress approaches its 
global failure point ac from below: x ~ ("^c —(^av)'^^"^ and r ~ {gc —cTav)"^^'^ ■ These results 
for the RFB model are of course analytically derived here for uniform distribution of strength 
of the fibers. It may be mentioned here that a similar behaviour for the time-to-fracture (for 
a above ac] diverging with the same exponent 1/2 for r) was observed in a RFB model where 
the fibers relax, under stress, to the elastic strain through viscous damping [2^. However, 
the relaxational dynamics in this visco-elastic RFB model is not due to the (self-organising) 
stress redistributions among the surviving fibers and, as such, is quite different in its origin. 
In fact, this time-to-failure vanishes in the limit of zero damping coefficient [pO|| . However, 



the similarities in the behaviour in such distinctly diff'erent situations also indicate interesting 
possibilities. 

Knowledge of the precursors and their power laws should help estimating precisely the 
location of the global failure or critical point from the proper extrapolation of the above 
quantities, which are available long before the failure occurs. The usefulness of such precur- 
sors can hardly be overemphasized. 
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Figure captions 

Fig. 1: The growth of average height hav (< hc{L)) of the BTW model against the number 
of iterations of adding unit heights (L = 100) . In the inset, we show the finite size behaviour 
of the critical height hc{L), obtained from simulation results for different L. 
Fig. 2: The variations of the precursors with hav (< hc{L)) in the BTW model for different 
system sizes: L = 100 (plus) L = 200 (cross) and L = 300 (open circle), (a) For relaxation 
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time r; in the inset r~°'^ is plotted against hav (b) For the total number of topplings A; 
inset shows A~°-^ versus hav plot, (c) For the correlation length ^; in the inset, ^^^'^ is 
plotted against hav 

Fig. 3: The growth of average height hav (< hc{L)) of the Manna model against the 
number of iterations of adding unit heights (L = 100) . In the inset, we show the finite size 
dependence of the critical height hc{L), obtained from simulation results for different L. 
Fig. 4: The variations of the precursors with hav (< hc{L)) in the Manna model for different 
system sizes: L = 100 (plus) L = 200 (cross) and L = 300 (open circle), (a) For relaxation 
time r; in the inset r~°-^ is plotted against hav (b) For the total number of topplings A; 
inset shows A~°-^ versus hav plot, (c) For the correlation length ^; in the inset, ^^^^ is 
plotted against hav 

Fig. 5: Fraction of the unbroken fibers Un at different times or iterations n in a RFB model 
with uniform strength distribution, for different values of (initial) stress: a = 0.24 (plus), 
a = 0.245 (cross), a = 0.248 (open circle), a = 0.25 (open square), a = 0.252 (open triangle). 
Note that the last value of a is greater than o-c{= 1/4), and the fraction of unbroken fibers 
goes to zero here. Inset shows how the susceptibility x (up triangle) and the relaxation time 
T (filled circle) both diverge as a — * o"c(= 1/4)- 
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Fig. 3 
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